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ABSTRACT 

We discuss non-self-gravitating hydro dynamic disks in the thin disk limit. 
These systems are stable according to the Rayleigh criterion, and yet there is 
some evidence that the dissipative and transport processes in these disks are hy- 
drodynamic in nature at least some of the time. We draw on recent work on the 
hydrodynamics of laboratory shear flows. Such flows are often experimentally un- 
stable even in the absence of a linear instability. The transition to turbulence in 
these systems, as well as the large linear transient amplification of initial distur- 
bances, may depend upon the non-self-adjoint nature of the differential operator 
that describes the dynamics of perturbations to the background state. We find 
that small initial perturbations can produce large growth in accretion disks in 
the shearing sheet approximation with shearing box boundary conditions, despite 
the absence of any linear instability. Furthermore, the differential operator that 
propagates initial conditions forward in time is asymptotically close (as a func- 
tion of Reynolds number) to possessing growing eigenmodes. The similarity to 
the dynamics of laboratory shear flows is suggestive that accretion disks might be 
hydrodynamically unstable despite the lack of any known instability mechanism. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities 
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1. Introduction 

In this paper we discuss the hydro dynamic stability of thin, non-self-gravitating, 
Keplerian accretion disks. The study of potential hydrodynamic and magnetohydrodynamic 
(MHD) instabilities in Keplerian shear flow is important because they may explain transport 
within accretion disks. (For a review of accretion disk physics, see Frank, King, & Raine 
(12). A more recent review of accretion disk stability and transport may be found in Balbus 
& Hawley (3).) 

The effective kinematic viscosity v of accretion disks is traditionally parametrized by a 
in the phenomenological relation 

v = ac s H, (1) 

where c s is the sound speed and H is the local disk scale height (26). Inferred values of 
a for dwarf nova (DN) accretion disks, which are relatively well-studied and for which a 
is perhaps the most observationally constrained, lie in the range of a 0.1 — 0.001. (See 
Warner (30) for an excellent review of DN disks.) These values for a are many orders 
of magnitude (e.g., 10-12) larger than back-of-the-envelope estimates for the molecular 
viscosity. Other disk systems, such as the accretion disks in active galactic nuclei (AGN) 
and the circumstellar disks in young stellar objects (YSO), are also thought to have effective 
viscosities that are anomalously large in the same sense (21; 17). 

It was originally hoped that the large Reynolds numbers in Keplerian disks would be 
a source of hydrodynamic turbulence (29). However, Keplerian disks satisfy the Rayleigh 
stability criterion for inviscid axisymmetric perturbations (5), 

4(r 2 fi) 2 > 0. (2) 
or 

It should be noted that this is a necessary, but not sufficient, criterion for stability of the 
more general case of non-axisymmetric perturbations to a viscous flow. Nonetheless, there 
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is no well-established instability mechanism in Keplerian hydrodynamic shear flow. 

The canonical route to turbulence in disks is the Balbus-Hawley magnetic shearing 
instability (2; 18). However, there may be some disks that are not sufficiently ionized for 
the Balbus-Hawley instability to be effective. Ion-neutral collisions in circumstellar disks 
associated with YSOs may dissipate the magnetic field too quickly for the instability to 
take hold (27). Gammie & Menou (13) simulated the dynamics of the prototypical DN 
system SS Cyg, using a code developed by Hameury et al. (16); based on their simulations, 
they conclude that the Balbus-Hawley instability may not be operating when the disk is in 
quiescence. In their disk simulations, the magnetic Reynolds number drops below 10 4 , and 
numerical MHD simulations (19) indicate that the Balbus-Hawley mechanism requires a 
magnetic Reynolds number Res based upon scale height to be larger than approximately 
10 4 . 

We show that the shearing-sheet approximation with shearing-box boundary conditions 
shares properties in common with such laboratory flows as Couette flow, Poiseuille flow, 
Hagen-Poisseuille flow, and Taylor-Couette flow. (For an introduction to these flows, 
see Drazin & Reid (8).) These shear-dominated laboratory flows exhibit transitions to 
turbulence even in the absence of a linear instability (28; 24). Hagen- Poiseuille flow (pipe 
flow), for example, does not possess a linear instability, but is seen to become turbulent in 
the laboratory at a range of Reynolds numbers, as observed by Reynolds (23). 

One property of these flows (10; 28; 14; 22; 4; 15) is that initial perturbations may 
undergo large transient growth in the absence of any growing eigenmode. Furthermore, the 
eigenvalue spectrum of these operators is very sensitive to perturbations to the operators, 
so that the operators are in many cases asymptotically close (as a function of Reynolds 
numbers) to operators that possess growing eigenmodes, even if the operators themselves 
do not have growing modes. 
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2. Local hydrodynamics and shearing box 

The shearing-sheet approximation replaces the (r, 6) coordinates in the disk with a 
local Cartesian plane (x, y) (see e.g. page 439 of Ryu & Goodman (25) and references 
therein). We neglect the z dimension and assume incompressibility. We have then 

d t v + v • Vv + p- 1 Vp + 2fiz x v - AAVtx± = z/V 2 v, (3) 

where A is Oort's first galactic constant, 



-\rd r Q 



= f flo- (4) 

r 



The background flow is — 2Axy. Perturbations u to this flow may be described in terms of 
a streamfunction tp, so that v = — 2Axy + u, and (u x ,u y ) = (d y ip, —d x ip). We let 

^(x,y,t)=((>(x,t)e i ^, (5) 

and we non-dimensionalize by measuring time in units of 1/Q and distance in units of some 
unspecified L, which we may take to represent, say, the disk scale height. The dimensionful 
viscosity v is the inverse Reynolds number v — > R^ 1 with the length L as the defining 
lengthscale. In the case that L is the disk scale height, then R may typically lie in the range 
10 10 -10 12 approximately, as discussed above. 

Taking the curl of (3) and neglecting d z u z in the spirit of a 2D analysis, we have then 

that 

d t 4> = [dl - [3 2 ] - 1 ip [(iPR)-\% - P 2 ) 2 - 2Ax(d 2 x - P 2 )] 0, (6) 

which is a form of the Orr-Sommerfeld equation for viscous shear flow (Reddy et al. 1993). 
(Note that in Keplerian disks, the nondimensionalized A is simply A = 3/4.) If we further 
separate variables, this equation is seen not to be formally self-adjoint in x in the L 2 norm 
for 0. With the addition of rigid boundaries in x, this equation governs 2D perturbations 
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to Couette flow; at high Reynolds number, its eigenmodes are highly collinear, and its 
eigenvalues are extremely sensitive to perturbations to the equation. 

We seek solutions to equation (6) in the shearing-box system. This is the shearing 
sheet with the addition of periodic boundary conditions in y and sheared periodic boundary 
conditions in x (3). Specifically, we consider a square domain of size L, and the boundary 
condition in x is 

if>(0,y,t) = il>(L,y-2ALt,t). (7) 

Equation (6) is not separable in the shearing box system. However, a complete set of 
solutions exist; these are the Kelvin modes (20; 11). These modes are Fourier modes in x 
and y, where 

^(x,y,t) = C m ^(t)e ik ^ +i ^, (8) 

with k{t) = fc(0) + 2 Apt, and 

M0) 2 + /3 2 „„„ o-i / ,ur^2 . ^ (g) 



Ck^At) = C H0)A°) UU \2 i Q2 eXP 



-R- 1 I (k(r) 2 + (3 2 ) dr 
Jo 



k(t) 2 + P 

To satisfy boundary conditions, we have (5 = m2n and k(0) = n2n. At times t = jAt, 
where j is an integer and At = (2Am) _1 , the streamfunction may be written 

^(x,y,t) = cW m e 2 ^ nx+m «\ (10) 

We restrict our attention to the invariant subspace m — 1, since these are the least 
dissipated modes with the potential for growth, and we henceforth drop the subscript m. 
The streamfunction at timestep j is then completely specified in terms of the state vector 
c^'l = (..., c^, Cq \ c± \ . . .). The solution to equation (6) may now be written as a difference 
equation in t, 

c \j+i]=CcW, (11) 

where £ is a matrix operator. It is zero everywhere except on the subdiagonal, reflecting 
the fact that the i th element of the state vector is taken from the (i — l) th element of 
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c^ 1 ' (i.e. c on the previous timestep), adjusted in amplitude according to equation (9). 
The matrix £ provides a convenient representation of the propagator for the initial value 
problem: 

C = £ j S (12) 

The operator £ is formally generated by the Orr-Sommerfeld operator £qs (which is the 
right-hand side of equation 6), in the subspace of discrete Kelvin modes described above: 

£ = e AtCos . (13) 

The Laplacian in the operator £ Q s is rendered well-defined by the boundary conditions. 



3. Transient growth and operator perturbations 

In order to discuss growth, we first need a norm to measure the strength of 
perturbations. We use the L 2 norm of the streamfunction ip. In this norm, the magnitude 
of a perturbation c is 




if)* if) dx dy = c*q. (14) 



Given this vector norm, the usual operator norm of a matrix A is defined using the vector 
norm: 

= max| N | =1 ||Ar||, (15) 

E1 and it is equal to the maximum singular value of A. The norm measures the 

maximum possible growth of perturbations in j timesteps. 

In figure 1, we plot (triangles) the operator norm ||>C'|| versus j, for a Reynolds number 
of 10 6 . The "hump" is a typical characteristic behavior for the norm for a non-normal 



E1 NOTE TO EDITOR: The output for equation 15 is not quite what I would like: The 
"max" needs to be directly over the "||x|| = 1", but I don't know how to do that. 
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matrix. We show in figure 2 (triangles) the log of the amount of growth that can be obtained 

as a function of Reynolds number, which we denote by G(R), where G(R) = max 

and the maximum is taken over j. As C is represented by a matrix that is zero everywhere 

except on the subdiagonal, powers of it are easy to compute. The maximum growth scales 

as 

G(R) ~ f? 2/3 , (16) 

where the maximum is taken over the values of j. This result is also easily obtained 
analytically from equation (9) in the limit of large R. Specifically, for large R, we have 

G(R) « (0.092)i? 2 / 3 (17) 

The line drawn through the triangles in figure 2 is the analytic approximation above. 

This type of transient growth has been remarked upon previously in the literature; it is 
merely a consequence of the kinematics of the Kelvin modes as described by equation (9). 
According to the Rayleigh hypothesis, instability occurs whenever there exists an unstable 
mode. The unbounded exponential growth of an unstable mode is of course unphysical. 
Large transient growth in the linear limit, as obtained here, simply implies that the strength 
of the initial perturbation required in order to reach the nonlinear regime is potentially very 
small, but not infinitesimal. Nonlinear effects may then positively or negatively feedback 
upon the original perturbation, but we do not discuss these effects here. 

This transient growth has implications for the spectra of perturbations to the operator 
C. We construct matrices M. such that each element of M. is randomly drawn from the 
unit disk in the complex plane. We then construct the perturbed operator C, 

C = C + eM. (18) 

The matrix M. couples the Kelvin modes, allowing the possibility for feedback to the 
transient growth observed in ||£ n ||. Note that we do not rescale M. by its operator norm, 
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||M|| « y^iMM)- 

Although in principle M. need not have any eigenvalues with magnitude larger than 
one, in practice it always will. Given this, for each M. there exits a positive e such that 
the operator C possesses an eigenvalue A with |A| > 1, so that the operator C possesses a 
growing eigenmode. Such an eigenmode corresponds to a growing Bloch mode, 

<t>(x,t)=t;(x,t)e i « t , (19) 

where t) is periodic in t with period At. In figure 1 (open squares), we show the growth 
of the operator norm | \(C + eA4) n \ \ at a Reynolds number of 10 6 and with e = 5 x 1CT 4 . The 
behavior of this perturbed operator diverges exponentially from that of the unperturbed 
operator in the operator norm, at large times. This reflects the presence of a growing 
eigenmode. Interestingly, the minimum values of e (denoted e crit ) required to produce a 
growing eigenmode for the operator C do not depend greatly upon the matrix A4, for our 
choice of probability measure for M.. Furthermore, the values e cr i t become asymptotically 
small as the Reynolds number is increased. 

In figure 2 (circles) we plot the log of 1/ e crit versus the log of the Reynolds number. 
We perform this calculation for forty different random matrices for each Reynolds number. 
(This is a more computationally intensive operation than the determination of the maximum 
operator norm maxU^'U. We have made extensive use of the LAPACK routine package 
(1), running on the beowulf cluster parallel virtual machine in the University of Texas 
Astronomy Department.) The points are the median values of 1 / e C rit, and the error bars are 
the first and third quartile. The relationship is a power law with an index of 0.88 ± 0.01. 
Deviations from this power law at the low end reflect the effects of the discretization of 
time into timesteps At. Deviations from the power law at high Reynolds number are due to 
computational limits on the necessarily finite size of the matrix representations for £ and 
A4, as we have confirmed by varying these dimensions. The corresponding points for larger 
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or smaller matrices fall on top of the points shown here, up until the point at which they 
turn off of the power law relationship. This point is determined by the value for j = j pca k at 
which the operator norm U^'H reaches its maximum, as matrices of order smaller than j will 
not faithfully represent the growth of ||£ J, ||. In figure 2 we also show the dependence of j pe ak 
upon Reynolds number (squares); it is seen to depend upon R in the form n ~ R 1 / 3 , and 
the line is the theoretical result for large Reynolds number. Given maximum growth G(R) 
of the operator norm for the unperturbed operator, we need only perturb the operator by 
a matrix that is zero everywhere except at the matrix element that couples Kelvin modes 
at the end of their growth period to the modes that are just beginning their growth period 
, i.e. jpeak timesteps earlier. This element need only have an amplitude of e > 1/G(R), to 
create a feedback loop. Rougly, the number of matrix elements that will effectively couple 
grown modes to the earlier maximally amplified modes may be expected to scale as jp eak ; 
when added in quadrature, the effective feedback then would scale as j pea kG(R), which 
grows as R 10 . However, each of these matrix elements do not have the same potential 
for producing feedback to the Kelvin modes, so that the actual scaling of G(R) with R is 
somewhat less steep. 

One physical interpretation of the operator perturbation eM. is that it potentially 
represents the differences between the idealized operator C and the true operator that 
governs the system. Compressibility, tides, magnetic fields, and so forth, will affect the 
operator C. While we cannot say that any of these effects will actually result in a new 
operator C that possesses growing modes, we can say that, at least in the probability 
measure of random matrices Ai that we have discussed, most perturbed operators of the 
form C = C + eM. have growing eigenmodes for asymptotically small e at high Reynolds 
number. This result reflects properties of C itself. Had C been a normal operator, this result 
would not obtain, and the spectrum of C would be very stable to operator perturbations of 
the form described above. 



4. Conclusions and Discussion 



The shearing sheet system with shearing box boundary conditions appears, as 
hypothesized, to share many properties in common with such flows as Couette flow and 
other shear-dominated flows. The operator C discussed above is clearly not normal. It does 
not even possess eigenmodes. However, as we have shown, C is capable of large transient 
growth in the operator norm, and it is within e of possessing growing eigenmodes in the 
sense described above, where e ~ Rr 1 . In particular, for the range of Reynolds numbers 
expected in disks as discussed above, e is on the order of 10~ 8 5 . 

This suggests that Keplerian shear flow may share the same phenomenology as the 
above-mentioned classic laboratory flows, including turbulence in the nominally stable 
regions of parameter space. We hypothesize that non-magnetized accretion disks may 
therefore be turbulent in the absence of a linear instability, as suggested previously by 
Dubrulle (9) and Richard & Zahn (24). This could explain the presence of transport in 
DN in quiescence, as well as transport in other cool disks such as those associated with 
young stars. We suggest that simulations of 2D vertically-symmetric hydrodynamics in 
the shearing box system should exhibit large transient growth of initial perturbations and 
operator sensitivity as described above. 

The author wishes to thank the Theory Group of the Department of Astronomy at the 
University of Texas for their continued support, and for their generosity with the group 
computing facilities. He also wishes to thank Ralph Showalter of the Math Department for 
many helpful discussions. 
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Fig. I. — II^H vs. j for R = 10 6 . Triangles are the unperturbed operator, squares are for 
the operator perturbed by a random matrix with matrix elements taken from the unit disk 
and scaled by e = 5 x 10~ 4 . Resolution is 201 modes. 

Fig. 2. — Log-log plot of e~; t (closed and open circles), G(R) (triangles), and j pca k (squares) 
versus Reynolds number R. Dashed line and dotted line are theoretical curves for large R 
limit for G(R) and j pea k respectively, and solid line is power law fit to the central part of 
curve for e~; t . Deviation from power law behavior at high R is due to the limitations of 
numerical resolution, as discussed in text. Closed and open circles are results from runs with 
401 and 801 modes respectively. 
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5. POSTSCRIPT 

The document appearing on the previous pages was submitted to ApJL on 19 Oct 2000 
and appears here now without alteration. This paper was not accepted for publication. 
The primary justifications for this decision were that (1) the streamfunction norm is not a 
physical norm, (2) the measure of perturbations in the form of random matrices as described 
does not have an obvious physical basis, (3) it was suggested that the 2D perturbations 
presented here are essentially uninteresting and that the perturbations should be performed 
in 3D. Note, for example, that the perturbations that experience the greatest transient 
amplification in planar Couette flow consist of shear-aligned vortices. These are inherently 
three-dimensional structures, Squire's theorem not being relevant. However, 3D streamwise 
vortices in Keplerian shear flow are quite different from similar perturbations to plane 
Couette flow, because of the presence of an epicyclic frequency, so to an extent this is where 
the interesting physics is. 

It was also pointed out that the bootstrapping described here "cannot be regarded as 
a proof of nonlinear instability." On that the referee and I are wholly agreed. If I implied 
in my paper that I thought differently, it was unintentional. 

In any case, the first objection listed above is quite valid, but on the other hand this 
problem is easy enough to fix, and switching to a more physical norm such as the energy 
norm has little material effect on the results. 

The second objection is also valid. The problem is that, as the referee put it, "a 
random perturbation of a system of equations embodying a conservation law will allow an 
artificial growth of the previously conserved quantity." True. (Although, strictly speaking, 
here the unperturbed system of equations is dissipative.) This defect might be remedied 
by restating the perturbations in terms of, say, noise in the form of random forcing in the 
system. Putting the perturbations on more physical grounds would have made the paper 
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much better, and probably would not have required too much work in the end. 

The third objection is quite obviously more substantial in the obstacles it presented to 
me. When I submitted this paper, I had already accepted employment elsewhere working 
in a completely different subject area. I felt I could fix the other problems raised by the 
referee, but a 3D analysis would take substantially more work; without question I did not 
have the time to do such an analysis in place of the 2D results presented here. 

I was at the time completely naive about the publication process. This was the first 
paper I had ever submitted, as is evident by my lack of familiarity with the idiom. Although 
I was very upset by the report, I simply deferred to the referee on this matter. I reluctantly 
abandoned the paper. Perhaps he was right, I thought, and only a 3D analysis would be 
worthy of publication in the peer-reviewed literature. 

I sincerely wish to thank the referee for the time and effort he took for his quite 
thorough and penetrating analysis. In reading the report over now, I see that he was trying 
to make helpful and encouraging suggestions, but all I heard at the time was a flat "no." 
Understanding the editorial process, and reading between the lines of referee reports, it 
appears, are acquired arts. On the referee's third objection above, however, I must now 
respectfully disagree. With perhaps a fair bit of polish, following the remaining suggestions 
of the referee, this modest work might have been a helpful addition to the literature at the 
time. I believe now, as I did in October of 2000, that even a 2D analysis of this system was 
interesting and worthy of appearing in the peer-reviewed literature. 

Apparently I am not alone. In the intervening years, other researchers have 
performed similar, although admittedly more insightful and extensive, analyses of transient 
amplification in this 2D system. I mention a few noteworthy papers below, some of which 
were accepted for publication by the referees and editors of other journals. My list of 
relevant papers is certainly not exhaustive and I apologize to anyone I have left out. I 
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have been quite happy to see these papers come out, because I believe in the value of this 
type of analysis. It has been illuminating to see the different approaches that people have 
taken to this problem. I anticipate that this line of work will ultimately lead to a greater 
understanding of transport in cool disks. 

There are many stories behind my paper, each with its own amusing ironies. These in 
turn mostly serve now as humorous little lessons to me that I will largely keep to myself. As 
I have no intention whatsoever of resuming the work presented here, I have decided simply 
to post it for archival purposes. A more extended discussion appears in the second chapter 
of my dissertation. So, if anybody wonders why I did not publish my dissertation, here at 
least is part of the answer. :) 

My graduate research advisors, J. Scalo and J. C. Wheeler, deserve my ample thanks 
for their support. My work certainly did not benefit either of their research programs; I 
had conceived of the project long before I ended up under their wings. In retrospect, it was 
clearly a mistake on my part to continue with a project I had begun earlier in my career, 
but "it seemed like a good idea at the time," as they say. As I worked autonomously on 
this project, any mistakes or other gaffes in this paper are entirely my own. 

Thanks go to Matt Umurhan for pointing out to me several of the references listed 
below, and for some fun and enlightening conversations on related topics. I also thank the 
reader for enduring my effusive commentary. I felt I couldn't post this paper without an 
explanation. 

Peter Williams 
Albany, CA 



-19- 
REFERENCES 

Ashfordi, N., Mukhopadhyay, B., & Narayan, R 2005, ApJ, 629, 373, arXiv:astro-ph/0412194 

Chagelishvili, G. D., Zahn, J.-R, Tevzadae, A. G., & Lominadze, J. G. 2003, A&A, 402, 
401, arXiv:astro-ph/0302258 

Mukhopadhyay, B., Ashfordi, N., & Narayan, R. 2005, Proceedings of the 22nd Texas 
Symposium on Relativistic Astrophysics, Dec 12-17, 2004, arXiv:astro-ph/0501468 

Umurhan, O. M., & Regev, O. 2004, A&A, 427,, 855, arXiv:astro-ph/0404020 

Yecko, P 2004, A&A, 425, 385. 



This manuscript was prepared with the AAS IATjtX macros v5.2. 



